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Abstract 

An assumed-stress hybrid/mixed 4-node quadrilateral shell element is introduced that al- 
leviates most of the deficiencies associated with such elements. The formulation of the 
element is based on the assumed-stress hybrid/mixed method using the Hellinger-Reissner 
variational principle. The membrane part of the element has 12 degrees of freedom in- 
cluding rotational or “drilling” degrees of freedom at the nodes. The bending part of 
the element also has 12 degrees of freedom. The bending part of the element uses the 
Reissner-Mindlin plate theory which takes into account the transverse shear contributions. 
The element formulation is derived from an 8-node isoparametric element. This process 
is accomplished by assuming quadratic variations for both in-plane and out-of-plane dis- 
placement fields and linear variations for both in-plane and out-of-plane rotation fields 
along the edges of the element. In addition, the degrees of freedom at midside nodes are 
approximated in terms of the degrees of freedom at corner nodes. During this process the 
rotational degrees of freedom at the corner nodes enter into the formulation of the element. 
The stress field is expressed in the element natural-coordinate system such that the element 
remains invariant with respect to node numbering. The membrane part of the element is 
based on a 9-parameter stress field, while the bending part of the element is based on a 
13-parameter stress field. The element passes the patch test, is nearly insensitive to mesh 
distortion, does not “lock”, possesses the desirable invariance properties, has no spurious 
modes, and for the majority of test cases used in this paper produces more accurate results 
than the other elements employed herein for comparison. 

Introduction 

The finite element method is a powerful technique that is used to solve a variety of very 
complicated problems in different fields of engineering and science. In the field of struc- 
tural mechanics, the finite element method is the most widely used method for obtaining 
solutions to structural analysis problems. The finite element modeling of general aerospace 
shell structures usually requires the use of distorted meshes. Additional complexities arise 
if the original finite element discretization is further refined through an adaptive refinement 
procedures. Most of the 4-node shell elements developed in the past do not produce reliable 
results for distorted meshes. The element developers are usually aware of the limitations 
and pitfalls of the elements. The users, on the other hand, may not be aware of all these 
limitations and make invalid uses of the these elements. It is desirable to develop simple 
3-node and 4-node shell elements that are free from these limitations and pitfalls, such as 
spurious zero-energy modes, locking, sensitivity to mesh distortion, invariant properties, 



and most importantly produce accurate and reliable results. The elements that have been 
developed to date lack one or more of these characteristics. The standard 3-node and 
4-node isoparametric elements are generally too stiff for most practical problems, and a 
“blizzard” of these elements must be used to model even simple shell structures in order 
to obtain accurate and reliable results. In order to remedy the situation, researchers have 
developed many non-standard elements which have their own merits and shortcomings. 
Some of these elements do not pass the patch test — a requirement for convergence to the 
correct solution. Some other elements have spurious mechanisms (zero-energy deforma- 
tional modes) such as the hour-glass mode, while some others do not possess the desirable 
invariance properties with respect to node numbering. Furthermore, most elements are 
sensitive to geometry or mesh distortion. 

One way to remedy the deficiencies of the membrane elements is to include the nodal ro- 
tational or “drilling” degrees of freedom in the element formulation. In early attempts, 
these rotational degrees of freedom were used in cubic displacement functions. However, 
Irons and Ahmad demonstrated that this concept had serious deficiencies [1]. The elements 
formed in such a manner force the shearing strain to be zero at the nodes and do not pass 
the patch test, which could produce serious error in some structural analysis problems. 
Recently several researchers have used these rotational degrees of freedom in quadratic 
displacement functions with more success[2-5]. This latter method is constructed in the 
following way. First, the element is internally assumed to be an 8-node (4 corner nodes 
and 4 midside nodes) isoparametric element with 16 degrees of freedom and the stiffness 
matrix associated with this “internal” element is calculated. Then, this stiffness matrix is 
condensed to that of a 4-node element with 12 degrees of freedom by associating the dis- 
placement degrees of freedom at the midside nodes with the displacement and rotational 
degrees of freedom at the corner nodes. MacNeal[4] has used this approach to develop 
4-node displacement-based membrane elements with selective reduced-order integration. 
Yunus et al.[5] have also used this method to develop assumed-stress hybrid/mixed mem- 
brane e lements. 

In this paper, this method is extended to the plate problem as well as the membrane 
problem to develop a 4-node assumed-stress hybrid/mixed quadrilateral shell element with 
24 degrees of freedom. The formulation is based on the Hellinger-Reissner variational 
principle. The element formulation is derived from an 8-node isoparametric shell element 
with 40 degrees of freedom by eliminating the degrees of freedom at the midside nodes 
at the expense of adding the rotational degrees of freedom to the corner nodes. This 
process is accomplished by assuming quadratic variations for both in-plane and out-of- 
plane boundary displacements and linear variations for both in-plane and out-of-plane 
boundary rotations along the edges of the element and evaluating these functions at the 
midside locations. These midside values are then used to condense the stiffness matrix of 
the 8-node element to that of the 4-node element. As a result, the rotational degrees of 
freedom at the corner nodes enter into the formulation. 

The membrane part of the element has 12 degrees of freedom, thre e o f which account 
for the rigid body modes. Therefore, a minimum of 9 independent stress parameters are 
needed to avoid rank deficiency for the membrane part. The bending part of the element 
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also has 12 degrees of freedom, three of which account for the rigid body modes. Therefore, 
a minimum of 9 independent stress parameters are also needed to avoid rank deficiency for 
the bending part. The membrane stress field is represented with 9 independent parameters. 
The bending stress field, on the other hand, is represented with 13 independent parameters 
which was derived in references [6-8] through an asymptotic power-series expansion of 
stresses in the thickness direction. This 13 parameter selection of stresses for the bending 
part is less sensitive to geometry or mesh distortion than a 9 parameter selection obtained 
from a degenerate solid model [9]. 

Evolution of Hybrid Method 

The classical assumed-stress hybrid formulation of Pian[10] is based on the principle of 
minimum complementary energy. The displacements are described only along the element 
boundary, and an equilibrating stress field is described over the domain of the element. 
It was later recognized that the same method may be derived from the Hellinger-Reissner 
principle[ll-13]. However, in the Hellinger-Reissner principle the stress field is not required 
to satisfy the equilibrium equations a priori, while the displacement field must be described 
over the domain of the element and not just along the boundaries. The stress field would 
then satisfy the equilibrium equations only in a variational sense. Therefore, the stress field 
may be described in the natural-coordinate system of the element which would make the 
element less sensitive to mesh distortion, and a proper selection of the stress field would 
make the element invariant with respect to node numbering. For example, if the node 
numbering is changed from 1-2-3-4 to 2-3-4-1, then the natural coordinate of the element 
would rotate by 90°. Therefore, a selection of stresses that is invariant under a 90° rotation 
would make the element invariant. Because of these desirable attributes, researchers have 
developed assumed-stress hybrid/mixed elements using the Hellinger-Reissner principle. 

In this study, the Hellinger-Reissner principle is used to develop a 4-node shell element, 
which includes the rotational degrees of freedom in the membrane part of the element for- 
mulation. The details of the assumed-stress hybrid/mixed formulation using the Hellinger- 
Reissner principle has been extensively discussed in the literature (e.g;, see references [11- 
13]), and only a brief outline is given herein for completeness. The variational functional 
is given as 

n = -- / <x T DcrdV’ + / er T C[u]dV- [ u T t 0 dS (1) 

2 Jv Jv Js t 

where D is the compliance matrix of the material, er is the stress array, u is the displacement 
array, to is the prescribed traction array, matrix £ is a differential operator which produces 
a strain field when operated on the displacement field, V is the domain of the element, 
and S t is the part of the boundary where t 0 is specified. The assumed-stress hybrid/mixed 
formulation is based on assuming a stress field in the interior of the element as 

cr = P/3 (2) 

and assuming the displacement field as 

u = Nq (3) 
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where the matrices P and N consist of the appropriate interpolating functions Tor stresses 
and displacements, respectively, and the vectors (3 and q are the unknown stress parameters 
and nodal displacements and rotations, respectively. 

The expressions for stresses in equation (2) and displacements in equation (3) are substi- 
tuted into the functional II of equation (1) and the variation of the functional with respect 
to the element internal unknowns (3 are set to zero. This stationary condition gives 


where 


(3 = H -1 Tq 

( 4 ) 

H = / P t DP dV 

( 5 ) 

Jv 


T = J P T £[N]<fF 

( 6 ) 


The substitution of the expression for (3 in equation (4) into the functional II of equation 
(1) and a subsequent variation of the functional with respect to the global unknowns q 
yields 

Kq = F (7) 


where the stiffness matrix K is given by 


K = T t H _1 T 


( 8 ) 


and the generalized force vector F by 

F= f N T t 0 dS 

Js t 


( 9 ) 




Element Formulation 

Midside Displacement Approximation 

As discussed earlier, the stiffness matrix for the present 4-node element is obtained from 
the stiffness matrix associated with the “internal” 8-node element shown in Figure 1. This 
process is accomplished by associating the degrees of freedom at midside nodes with the 
degrees of freedom at corner nodes by expressing the displacement components along each 
edge of the element as quadratic functions in terms of the nodal displacements and rotations 
as follows. If a typical edge of a quadrilateral is considered to be a 2-node beam (see Fig. 1, 
edge 1) with constant shear deformations then a displacement component perpendicular 
to the beam axis, for example the v ! component, may be described as a quadratic function 
as 

v = a 0 + ai£ + a 2 £ 2 (10) 

where, £ is a non-dimensional coordinate in the direction of the axis of the beam such that 
node 1 of the beam is at £ = — 1 and node 2 of the beam is at £ = +1. 

4 


. mil if 



The rotation about the z axis, 9 Z , is given by 


6. 



( 11 ) 


where 7 ' is the constant shear deformation of the beam in the x-y plane. The four unknowns 
a 0 , ai, a 2 in equation (10), and 7' in equation (11) may be solved for in terms of the two 
nodal displacements v\ and v 2 and the two nodal rotations 9 Z \ and 0 z2 . The displacement 
component v' is given by 

.’ = 5(1 - M + |(1 + f)»i - 1(1 - - *.i ) (12) 

where l is the length of the beam. The normal or “drilling” rotation is given by 

9 Z — -(1 — £)9 Z 1 + -(1 + £)9 Z 2 (13) 


and the shearing strain by 


y (14) 

The displacement parallel to the beam axis is described as a linear function in terms of 
the two nodal displacements as 

»• = l(i -««; + l(i +«k (is) 

These concepts may be readily extended to triangular and quadrilateral elements. There- 
fore, the in-plane displacement and rotation fields on a typical edge of a quadrilateral 
element (e.g., edge 1 of Fig. 1) expressed in the element reference x-y coordinate system 
may be described as 

u = ^(1 - + 0 U 2 + - jp-(l ~ £ 2 )(^2 “ 0*0 

V = ^(1 - £)tq + ^(1 + £) u 2 g^-(l - £ 2 ){9z 2 - &z\ ) (16) 

0z = \( l-0«.l+|(l + 0«.2 

where, Axj and Ayi are the Ax and A y of edge 1 with respect to the reference local 
element x-y coordinate system (e.g., Axj = x 2 — xj). Note that only the displacement 
functions described in equations (16) enter into the membrane formulation. It should be 
mentioned here that the true nodal rotations are given by ^(^r — §^) evaluated at the 
nodes. Hence, the terms 9 zi are not true nodal rotations, and they may be referred to as 
“rotational connectors” [ 2 ]. 
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In references [2-5], the above displacement functions were used for construction of mem- 
brane elements. However, this same procedure may also be used for development of a bend- 
ing element. Again, quadratic functions are used to describe the out-of-plane displacement 
component w and linear functions are used to described the out-of-plane rotations 0 X and 
0 y along each edge of the element. The final result for the out-of- plane displacement and 
rotation fields on edge 1 of Figure 1 is 

w = - 0 W 1 + + 0 W 2 - ^p(l - £ 2 )(0*2 - M + ^p-(l - £ 2 )(0y2 - Oyl) 


&x — -(1 — 0^*1 + 2(1 + 00x2 (17) 

+ 5(1 +0** 

Equations (16) and (17) indicate that both the membrane part and the bending part of 
the element are formulated in the same manner. All three displacement components are 
quadratic functions, while all three rotations are linear functions. This conformity in the 
order of the approximating polynomials for all displacement components and all rotation 
components is very desirable in the analysis of shell problems. 

The basic concept of the element is similar to that of Cook’s proposed 4-node membrane 
element[3]. However, this concept is extended in this paper to the bending problem in order 
to construct a 4-node flat shell element. Internally the element is an 8-node isoparametric 
quadrilateral element. Externally, the midside degrees of freedom of the 8-node element 
are approximated by the above displacement and rotation functions, so that the midside 
degrees of freedom are eliminated at the expense of adding drilling degrees of freedom to 
the 4- node element. For example, on edge 1 of Figure 1, the midside degrees of freedom 
are given by evaluating the functions in equations (16) and (17) at £ = 0. This gives 

u s — 2^ Ul U2 ) ^ S~^ z2 ~ 

V5 = ^(Vl + V 2 ) - -^-(0x2 ~ 0 Z 1 ) 

W 5 = + W 2 ) - ~~~{0 X 2 ~0 xl ) + ^p-(0y2 ~ 0 y i ) (18) 

Ox 5 = 2 (Oxl + 0 X 2 ) 

0y5 — ~{0y 1 + 0y 2 ) 

The description for other midside nodes is readily obtained. It is not necessary to evaluate 
the function 6 Z in equations (16) at £ = 0 since the normal rotation 0 Z does not enter 
into the membrane formulation. It is observed that the extra nodal degrees of freedom are 
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associated with the membrane part of the element only. The bending part of the problem 
does not add any extra degrees of freedom to the element. 

Equations (16), when extended to all four sides of the element, indicate that the membrane 
part of the element has 12 degrees of freedom. Two of these are the in-plane translational 
rigid body motions and one is the in-plane rotational rigid body motion. Of the nine 
remaining degrees of freedom, three represent the constant strain states, five represent 
higher-order strain states, and the final degree of freedom represents a special type of zero- 
energy “spurious” mode. This zero-energy mode is associated with a state of zero nodal 
displacements and equal nodal rotations. This state renders the u and v displacement 
components in equations (16) zero on all four sides of the quadrilateral, while all four 
corner nodes of the quadrilateral are rotated about its normal by equal amounts 0 Z . This 
mode is shown in Figure 2 using cubic interpolation for displacements and may be called 
a “zero displacement” mode[2]. It occurs because the quadratic displacement functions 
are based on the differences of the nodal rotations and not the nodal rotations themselves. 
One of these differences is dependent on the other three. Therefore, the membrane part 
of the element has, in fact, only 11 independent displacement modes but is expressed in 
terms of 12 degrees of freedom. Hence, one of the normal rotational degrees of freedom 
is extraneous and must be eliminated. This zero-energy mode is of a special type and is 
different from other spurious mechanisms such as the hour-glass mode. In reference [4], this 
zero-energy mode is removed by artificially adding a small energy penalty to the stiffness 
matrix to make the stiffness matrix non-singular. Although, this could be accomplished 
with little effort, prescribing the value of normal rotation for at least one node in the entire 
finite element model of the structure will eliminate this extraneous degree of freedom. 

Equations (17) indicate that the bending part of the element also has 12 degrees of freedom. 
One of these is the out-of-plane translational rigid body motion and two of these are the 
out-of-plane rotational rigid body motions. Of the nine remaining degrees of freedom, 
three represent the constant curvature states, two represent the constant transverse shear 
strain states, and the other four represent higher-order strain states. No zero-energy modes 
are associated with the bending part of the element. Note that although the displacement 
function w described in equations (17) is also in terms of the differences of the nodal 
rotations and not the rotations themselves, each one of the nodal rotations are accounted 
for by the functions 9 X and 9 y in equations (17). This is not the case for the membrane 
part because the rotation function 9 Z in equations (16) does not enter into the membrane 
formulation, while both rotation functions 9 X and 9 y in equations (17) enter into the 
bending formulation. 

Stress Field Description 

The stress field should be selected in such a manner that no spurious zero-energy mode 
is produced. A spurious zero-energy mode is produced when the product of a selected 
stress term and the strains that are derived from the displacement functions produces 
zero strain energy under a particular, but not trivial, deformational displacement field. In 
order to avoid spurious zero-energy modes, each independent stress term must suppress 
one independent deformation mode. The problem of spurious zero-energy modes generally 
occurs for regular geometries such as rectangular planar elements and brick solid elements 
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and it disappears for irregular geometries[14]. 

As discussed previously, the membrane part of the element has 12 degrees of freedom — 
three of which are due to the in-plane rigid body motions. Therefore, a stress field with 
a minimum of 9 independent parameters is needed to describe the membrane stress (re- 
sultant) field. The following stress (resultant) field is first considered for the membrane 
part 

= 01 + 04V + 0e£ + 08£v 

N v = (3 2 + 0s£ + 07 v + (3g£i 7 ( 19 ) 

N< v =(h- 06V ~ 0i£ - \hv 2 - 

This stress (resultant) field is obtained by integrating, through the thickness, the com- 
ponents of the contravariant stress tensor expressed in the natural coordinates. These 
components must be transformed to the physical stress components in the element refer- 
ence coordinates by contravariant tensor transformation laws. However, the components 
of the transformation are evaluated at the origin of the natural coordinates in order to 
retain the constant terms in the stress field description[ll-13]. When the stress (resultant) 
field in equations (19) is specialized to the rectangular coordinate system, the stress resul- 
tants of the stress field assumed by Robinson[15,16] for his 8-node and 4-node membrane 
elements are obtained. In the above stress (resultant) field, the first five terms represent 
the stress field that was used in Pian’s 4-node assumed-stress hybrid membrane element 
which had no rotational degrees of freedom[ll-13]. The remaining terms are present to 
suppress the rotational degrees of freedom present in this formulation. However, experi- 
mentations with a single rectangular element revealed that the stress (resultant) field in 
equations (19) produces one spurious zero-energy mode. It should be mentioned here that 
Robinson’s suggested stress field worked for his formulation which was based on a cubic 
displacement field. Therefore, it is proposed herein to use the following 9-term stress field 
for the membrane part 


N( = 0i + 04 V + 06 £ + 06V 2 

Nr, = 02 + 0U + 0rv + 09? ( 2 °) 

Nfr =03 - 06V - 0r£ 

This stress field produces no spurious zero-energy mode for the assumed displacement field 
described in equations (16) for the membrane part. When the elements are rectangular in 
shape this selection of stress field satisfies the equations of equilibrium a priori. In reference 
[5], a different 9-term stress field was selected which assumes linear variation for all stress 
components which satisfies equilibrium in a variational sense for all element geometries 
including rectangular-shaped elements. 

The bending part of the element has 12 degrees of freedom — three of which are due to the 
out-of-plane rigid body motions. Therefore, a stress field with a minimum of 9 independent 
parameters is needed to describe the stress field. The following stress (resultant) field is 
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chosen for the bending part 


= (3\ -f 047] + /3g£ + 0$CV 

M v =p 2 + Pst + fav + 09tv (21) 

— 03 + p 10 £ + PllT] + 2^ 12 £ 2 + 2^ 13 ^ 2 

The transverse shear forces are taken to be derived from the moments such that the 
equilibrium equations are satisfied for rectangular elements 


Qt — M U + 

Qv = + M ViV 


which gives 

Qt = (P 't + 0u) + (fa + Pis) 7 ! ^ 3 ) 

Qv = (07 + 0io) + (09 + fil2)£ 

Again, this stress (resultant) field is obtained by integrating, through the thickness, the 
components of the contravariant stress tensor expressed in the natural coordinates and 
must be transformed to the physical components in the element reference coordinates by 
contravariant tensor transformation laws. As for the membrane part, the components of the 
transformation are evaluated at the origin of the natural coordinates in order to retain the 
constant terms in the stress field description. This stress field selection for the bending part 
satisfies the equations of equilibrium for rectangular-shaped elements a priori. When the 
above bending stress resultant field is specialized to the rectangular coordinate system, one 
obtains the resultants of the stress field that was derived in references [6-8]. The derivation 
of the stress field in references [6-8] consists of expressing the stress components as power 
series in the plate thickness, substituting these stresses into the continuum equations of 
elasticity, and equating the coefficients of like powers of the plate thickness. As discussed 
earlier, this 13 parameter selection of stresses for bending is less sensitive to geometric 
distortion than a 9 parameter selection obtained from a degenerate solid model[9]. This 
selection of stresses produces no spurious zero-energy modes. 

In this paper, the Reissner-Mindlin Plate theory is used for the bending part. The bending 
part is of class C° and takes into account the effects of transverse shear deformations by 
assuming constant transverse shear strains through the thickness of the plate. This means 
that the transverse shear stresses are also constant through the thickness of the plate. 
However, generally the transverse shear stresses are zero on the plate surfaces. Therefore, 
a parabolic variation of transverse shear stresses and strains through the plate thickness 
is more reasonable. To account for this discrepancy, a static correction factor of 5/6 is 
included in the transverse shear strain energy, see Reissner[17]. 

Other Element Matrices 

Other element matrices needed to complete the element formulation may be derived in a 
similar fashion. The element force vector F in equation (9) may be derived in a manner 
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consistent with the stiffness matrix. The surface traction and pressure loads are assumed to 
vary bilinearly over the element surface and the force vector F of equation (9) is calculated 
for the “internal” 8-node element. This force vector is then condensed to that of the 4-node 
element using the approximations for the midside degrees of freedom described earlier. For 
line load calculations, on the other hand, it is not necessary to use the “internal” 8-node 
element. The force vector F corresponding to line loads may be calculated directly based 
on the 4-node element by assuming linear variation of line loads on the edges of the element 
and using the displacement and rotation shape functions derived earlier. 

The mass and geometric stiffness matrices for the present element may also be derived 
from the “internal” 8-node element in the same manner as for the stiffness matrix. 


Numerical Results 

Assessment of the 4-node quadrilateral shell element developed in this paper is presented in 
this section. The element has been implemented in the CSM Testbed Software System[18] 
using the generic element processor template[19]. Selected test problems include the patch 
test, the straight cantilever beam, the curved cantilever beam, the twisted cantilever beam, 
the Cook’s tapered and swept panel, the thick-walled cylinder, and the Scordelis-Lo roof 
problem. The assumed-stress hybrid/mixed flat shell quadrilateral element derived in this 
paper will be referred to as AQR8 (Assumed-stress Quadrilateral Reduced from an 8-node 
element) in the following discussion for convenience. The results for the present element 
are compared with the results using the QUAD4 element of the MSC/NASTRAN from 
reference [20], the Q4S element of reference [4], the AQ element of reference [5], and the 
ES1/EX47, ES5/E410, and ES4/EX43 elements of the NASA Langley CSM Testbed. A 
brief description of these elements is given in the appendix. The dimensions and properties 
for the test problems are chosen in consistent units. 

Patch Test 

As the first test of the accuracy of the element, the patch test problem suggested in 
reference [20] is solved. This patch test is shown in Figure 3. Elements of arbitrary shapes 
are patched together to form a rectangular exterior boundary which makes it easy to apply 
boundary conditions corresponding to constant membrane strains and constant bending 
curvatures. The applied displacement boundary conditions and the theoretical solutions 
are also shown in Figure 3. The ability of the element to reproduce constant states of 
strains is an essential requirement for achieving convergence to the correct solution as the 
finite element mesh is refined. This requirement is observed by considering an individual 
element within a mesh with a complicated stress field. As the mesh is refined, the stresses 
within the elements tend towards a uniform value. Therefore, the elements that cannot 
produce a state of constant strains should not be trusted to converge to the correct solution 
as the mesh is refined[l]. The present element (AQR8) passes both the membrane and the 
bending patch tests with no error. The recovered strains and stresses are both exact. 

Straight Cantilever Beam 

As a second test, the straight cantilever beam problem suggested in reference [20] is solved 
for the three discretizations (six elements) shown in Figure 4. The constant and linearly 
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varying strains and curvatures are evoked by applying loads at the free end of the beam 
to test the ability of the element to recover these states of deformations. The theoretical 
results for extension, in-plane shear, out-of-plane shear, and in-plane moment are simply 
calculated from the elementary beam theory including shear deformations. The theoretical 
result for the twist is .03406, according to Timoshenko and Goodier’s Theory of Elastic- 
ity[21]. Reference [20] quotes the answer to be .03208. Analysis with three successively 
refined meshes converged to .03385 which is much closer to the theory of elasticity solu- 
tion than to that of reference [20]. Therefore, the solution from the theory of elasticity 
is taken herein for normalization purposes. The normalized results are shown in Table 1. 
The results show that all elements perform well for the rectangular mesh. However, for the 
trapezoidal and parallelogram meshes which contain considerable distortions, only the Q4S 
element and the present element (AQR8) are adequate. For the present element, while the 
error for all meshes and loads is less than 3.5%, the parallelogram mesh with a twist end 
load exhibits an error of 15.9% which shows the sensitivity of the selected bending stress 
field to mesh distortion. 

Curved Cantilever Beam 

Next, the curved cantilever beam problem shown in Figure 5 is solved. The beam is 
formed by a 90° circular arc. In-plane and out-of-plane loads are applied at the free end 
to produce in-plane and out-of-plane deformations, respectively. The theoretical solutions 
are taken to be those quoted in reference [20]. The normalized results from the present 
element (AQR8) are tabulated in Table 2. Also shown in Table 2 are the results from 
other elements for comparison. In this problem the mesh is distorted only slightly and 
the results for all elements are good. However, it is seen that the AQR8 element performs 
better than the other elements in the table, particularly for the in-plane shear loading. 

Twisted Cantilever Beam 

The twisted cantilever beam problem is shown in Figure 6. The beam is twisted 90° over 
its length. There are 12 elements along the length of the beam. Therefore, each element 
has a 7.5° warp. The beam is subjected to in-plane or out-of-plane unit loads at its free 
end. Both the membrane and the bending contributions are significant for either type 
of load. The theoretical solutions are taken to be those quoted in reference [20]. The 
normalized results from the present element (AQR8) are tabulated in Table 3. Also shown 
in Table 3 are the results from other elements for comparison. In this problem the mesh 
is comprised of rectangular shaped elements and there is no in-plane mesh distortion. 
Therefore, this problem tests the sensitivity of the elements due to warp distortion. The 
results show that, although the AQR8 element exhibits good performance, the QUAD4 
element performs even better for this problem. 

Cook’s Tapered and Swept Panel 

The tapered and swept panel with one edge clamped and the other edge loaded by a 
distributed shear force is analyzed next (see Fig. 7). This problem was used by Cook 
and many other researchers to test the sensitivities of finite elements due to geometric 
distortions. The panel was analyzed by a coarse 2x2 mesh and a finer 4x4 mesh. The 
reference solution for the vertical displacement at point C is taken to be 23.90 quoted in 
reference [5] as the best known answer. The normalized results for the present element 
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(AQR8) along with the results for other elements are shown in Table 4. The mesh for this 
problem is distorted only slightly and all the elements produce reasonable results. The 
AQR8 element, however, yields results that are closest to the best known answers. 

Thick- Walled Cylinder 

The thick-walled cylinder problem is shown in Figure 8. This problem was suggested by 
reference [20] to test the performance of finite elements for nearly incompressible materi- 
als. Plane strain conditions and radial symmetry are assumed so that the only non-zero 
displacement component is in the radial direction. The theoretical solution for this prob- 
lem is given in any standard text book on elasticity (e.g., see reference [21]) and is shown 
in Figure 8. The normalized results for the present element (AQR8) along with the re- 
sults for other elements are shown in Table 5. The mesh in this problem is distorted 
only slightly and the main purpose is to test the behavior of the elements for nearly in- 
compressible materials. It is observed that, while the performance of the hybrid elements 
(i.e., the ES4/EX43 and the AQR8 element) does not deteriorate as the Poisson ratio 
approaches 0.5, the other elements fail to retain their performance. This insensitivity to 
nearly incompressible materials is a trait of assumed-stress hybrid elements. 

Scordelis-Lo Roof 

Finally, the Scordelis-Lo roof shown in Figure 9 is analyzed. This structure is a singly 
curved shell problem in which both the membrane and the bending contributions to the 
deformation are significant. The result reported in most references is the vertical dis- 
placement at the midpoint of the free-edge. The theoretical value for this displacement is 
quoted in reference [22] to be 0.3086, but the reference value quoted in reference [20] is 
0.3024 for normalization of results. The latter value is also used herein for normalization 
purposes. Because of symmetry, only one quadrant of the problem is modeled. The mesh 
on one quadrant is chosen to be NxN for N=2,4,6,8,10 (N=number of elements along each 
edge) to show the convergence of the solutions for the AQR8 element. The results of the 
normalized displacement at the midside of the free-edge are shown in Table 6. The results 
for other elements are also shown in Table 6 for comparison. The mesh for this problem 
is comprised of rectangular shaped elements and all the elements in the table seem to 
perform well. 


Conclusions 

A 4-node quadrilateral shell element with 24 degrees of freedom has been developed which 
alleviates most of the deficiencies associated with 4-node shell elements. The element is 
based on the assumed-stress hybrid/mixed formulation using the Hellinger-Reissner princi- 
ple. The membrane part of this element has 12 degrees of freedom and includes the drilling 
(in-plane rotational) degrees of freedom at the nodes. The bending part of this element 
also has 12 degrees of freedom. The bending part is of class C° and takes into account the 
transverse shear deformations. Both in-plane and out-of-plane displacements are assumed 
to have quadratic variations along the edges of the element, while both in-plane and out- 
of-plane rotations are assumed to vary linearly. A 9-parameter stress field is assumed for 
the membrane part and a 13-parameter stress field is assumed for the bending part. The 
element formulation is derived from an 8-node isoparametric element by eliminating the 
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midside degrees of freedom in favor of rotational degrees of freedom at corner nodes. Al- 
though the concepts of the present element are simple, the derivation and implementation 
of the element is somewhat awkward and indirect. Nearly all element matrices are derived 
based on an “internal” 8-node element. 

Results from the numerical studies demonstrate that this element is accurate, passes both 
the membrane and the bending patch tests, does not “lock”, has no spurious modes, 
and is nearly insensitive to mesh distortion. The present element also has the desirable 
property of being invariant with respect to node numbering. The results indicate that 
whereas the present element offers little advantage for problems with nearly rectangular 
meshes, it produces more accurate results when dealing with distorted meshes, than the 
other elements for the majority of test cases in this paper. This is important in general 
modeling of shell problems and adaptive refinement. The results also indicate that the 
present element produces more accurate results than the other elements in this paper 
when dealing with nearly incompressible materials, which is a consequence of assumed- 
stress hybrid element formulation. 

Results obtained to date warrant further research in making the element formulation more 
direct and extending the capabilities of the element to include stability analysis, dynamic 
analysis, and non-linear analysis. 
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Appendix 

The following is a brief description of the elements used in Tables 1-6 for comparison with 
the element developed in this paper. 

The QUAD4 MSC/NASTRAN element is a 4-node isoparametric shell element with selec- 
tive reduced-order integration. The transverse shear uses a string-net approximation and 
augmented shear flexibility [23], This element was developed by MacNeal and is available 
in the MSC/NASTRAN finite element code. 

The Q4S element is a 4-node shell element in which the membrane part is formulated 
internally as an 8-node isoparametric element with selective reduced-order integration and 
later reduced to a 4-node element by eliminating the midside degrees of freedom in favor 
of rotational degrees of freedom at the corner nodes. This element was developed by 
MacNEAL[4], The bending part of the Q4S is the same as that of the QUAD4[4]. 

The AQ element is a 4-node assumed-stress hybrid/mixed membrane element which is 
formulated internally as an 8-node isoparametric membrane element and later reduced to 
a 4-node membrane element by eliminating the midside degrees of freedom in favor of 
rotational degrees of freedom at the corner nodes. This element was developed by Yunus 
et al.[5]. The only result reported in reference [5] for the cantilever beam problem in 
Table 1 using the AQ element is for the mesh with trapezoidal-shaped elements with a 
unit in-plane end moment. This result is reported to be .85. The difference in the results 
between the AQ membrane element and the membrane part of AQR8 is in the selection 
of the assumed-stress functions. 

The ES1/EX47 element is a 4-node C° isoparametric assumed natural-coordinate strain 
(ANS) shell element developed by Park and Stanley [24-25] and implemented in the CSM 
Testbed Software System[18] by Stanley using the generic element processor template[19]. 
This element is not invariant and does not pass the patch test. This element does not 
include the drilling degrees of freedom in the formulation. 

The ES5/E410 element is a 4-node C 1 shell element which was originally implemented in 
the STAGS finite element code and later in the CSM Testbed by Rankin[26], This element 
includes the rotational degrees of freedom in its formulation and uses cubic interpolation 
for all displacement fields. This element is not invariant and does not pass the patch test. 

The ES4/EX43 element is a simple 4-node C° isoparametric assumed-stress hybrid/mixed 
shell element implemented in the CSM Testbed by this author [27]. This element passes 
the patch test and is invariant with respect to node numbering. This element does not 
include the drilling degrees of freedom in its formulation and uses linear interpolation for 
all displacement and rotation fields. 
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Table 1. Normalized tip displacements in direction of 
loads for straight cantilever beam. 


Tip Loading 
Direction 

QUAD4 

MSC/ 

NASTRAN 

ESI/ 

EX47t 

ES5/ 

E410t 

ES4/ 

EX43 

AQR8 

(present) 

( 

a) rectangular shaped elements 

Extension. 

.995 

.995 

.994 

.996 

.998 

In-plane Shear 

.904* 

.904 

.915 

.993 

.993 

Out-of-Plane Shear 

.986 

.980 

.986 

.981 

.981 

Twist 

.941** 

.856 

.680 

1.023 

1.011 

End Moment 

— 

.910 

.914 

1.000 

1.000 

( 

]b) trapezoidal shaped elements 

Extension 

.996 

.761 

.991 

.999 

.998 

In-plane Shear 

.071* 

.305 

.813 

.052 

.986 

Out-of-Plane Shear 

.968 

.763 

# 

.075 

.965 

Twist 

.951** 

.843 

# 

1.034 

1.029 

End Moment 

— 

.505 

.822 

.102 

.996 

(c) parallelogram shaped 

elements 

Extension 

.996 

.966 

.989 

.999 

.998 

In-plane Shear 

.080* 

.324 

.794 

.632 

.977 

Out-of-Plane Shear 

.977 

.939 

.991 

.634 

.980 

Twist 

.945** 

.798 

.677 

1.166 

1.159 

End Moment 

— 

.315 

.806 

.781 

.989 


f These elements are not invariant and do not pass the patch test. 

* The results from MacNeal’s Q4S element for in-plane shear load are reported in ref- 
erence [4] to be .993, .988, and .986 for the meshes in (a), (b), and (c) respectively. 

** These results for twist were normalized with .03028 in reference |^20]. Herein, all 
the other results for twist are normalized using .03046 according to Timoshenko and 
Goodier’s Theory of Elasticity [21]. 

# The element produces a singular stiffness matrix for this mesh. 
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Table 2. Normalized tip displacements in direction of 
loads for curved cantilever beam. 


Tip Loading 
Direction 

QUAD4 

MSC/ 

NASTRAN 

ESI/ 

EX47 

ES5/ 

E410 

ES4/ 

EX43 

AQR8 

(present) 

In-plane Shear 

.833 

.929 

.938 

.888 

.997 

Out-of-Plane Shear 

.951 

.935 

.887 

.925 

.956 


Table 3. Normalized tip displacements in direction of 
loads for twisted cantilever beam. 


Tip Loading 
Direction 

QUAD4 

MSC/ 

NASTRAN 

ESI/ 

EX47 

ES5/ 

E410 

ES4/ 

EX43 

AQR8 

(present) 

In-plane Shear 

.993 

1.357 

1.054 

1.361 

.991 

Out-of-Plane Shear 

.985 

1.293 

1.173 

1.359 

1.093 


Table 4. Normalized vertical deflection at point C 
for tapered and swept panel. 


Mesh 

AQ 

ESI/ 

ES5/ 

ES4/ 

AQR8 



EX47 

E410 

EX43 

(present) 

2x2 

.914 

.880 

.873 

.882 

.930 

4x4 

.973 

.953 

.953 

.962 

.979 


Table 5. Normalized radial displacements at inner 
boundary for thick-walled cylinder. 


Poisson’s 

Ratio 

QUAD4 

MSC/ 

NASTRAN 


ES5/ 

E410 

ES4/ 

EX43 

AQR8 

(present) 

0.49 

.846 

.831 

.848 

.990 

.989 


.359 

.352 

.360 

.990 

.988 

0.4999 

.053 

.052 

.053 

.990 

.988 
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Table 0. Normalized displacements at the midpoint of 
the free-edge for Scordelis-Lo roof. 


Mesh 

QUAD4 

MSC/ 

NASTRAN 

ESI/ 

EX47 
_ 

m 

ES4/ 

EX43 

AQR8 

(present) 

2x2 

1.376 

1.387 

1.384 

1.459 

1.218 

4x4 

1.050 


1.049 


1.021 

6x6 

1.018 


1.015 


1.006 

8x8 

1.008 


1.005 


1.003 

10x10 

1.004 



1.011 

1.001 
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Location of nodes: 


node 

X 

y 

i 

.04 

.02 

2 

.18 

.03 

3 

.16 

.08 

4 

,08 

.08 


Applied displacements: 

(a) Membrane patch test 

Boundary conditions: u = 10 _3 (x + y/ 2) 

v = 10 _3 (x/2 +y) 

Theoretical solution: t xx — e yy = 7 X j/ = 10 3 

<7 XX — &yy = 1333., <T xy = 400. 

(b) Bending patch test 

Boundary conditions: w = — 10 _3 (x 2 + xy + y 2 )/ 2 

e x = — 10 _3 (x/2 +y) 

9 y = — 10 -3 (x + y/2) 

Theoretical solution: 

Bending moments per unit length: 

M x = M y = l.lllxlO' 7 , M xy = 3.333x10-® 
Surface stresses: 

a xx = <r yy = ±0.667, <r xy = ±0.200 

Figure 3. The patch test problem. a=0.24, b=0.12, t=0.001, E=10 6 , zv=0.25. 

(Consistent units are used for various properties.) 
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a) rectangular shape elements 


X 


A 


A 


45 


\ / A 


b) trapezoidal shape elements 


45 


/ / 


z^z 






c) parallelogram shape elements 


Theoretical solutions: 


Tip load direction 

Displacement in direction of load 

Extension 

.3xl0 -4 

In- plane shear 

.1081 

Out-of-plane shear 

.4321 

Twist 

.03406 

In-plane moment 

.009 


Figure 4. Straight cantilever beam problem. Length=6., height=0.2, depth=0.1, E=10 7 , 
v=0.3, mesh=6xl. Loading: unit forces at the free end. 
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Theoretical solutions: 

Tip load direction 

Displacement in direction of load 

In-plane shear 

.08734 

Out-of-plane shear 

.5022 


Figure 5. The curved cantilever beam problem. Inner radius=4.12, outer radius=4.32, 
depth=0.1, E=10 7 , i/=0.25, mesh=6xl. Loading: unit forces at the free end. 
(Consistent units are used for various properties.) 
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Tip load direction 

Displacement in direction of load 

In-plane shear 

.005424 

Out-of-plane shear 

.001754 


Figure 6. The twisted cantilever beam problem. Length=12., width=l.l, depth=0.32, 
twist— 90° (root to tip), E=29.xl0 6 , i^=0.22, mesh=12x2. Loading: unit 
forces at the free end. 

(Consistent units are used for various properties.) 
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Figure 7. The tapered and swept panel problem. Thickness=l., E— 1., v— 1/3, 
mesh=NxN. Loading: unit in-plane shear force distributed on the free edge. 
Reference solution: vertical displacement at C=23.90 from reference [10]. 
(a) 2x2 mesh, (b) 4x4 mesh. 

(Consistent units are used for various properties.) 
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Theoretical solutions: 


u(r) = 


(1 +v)pR\ 
E{R\ - R\) 


[Rl/r + ( l-2u)r] 


where p=pressure on the inner surface, Ri dinner radius, i 22 =outer radius. 


Poisson’s ratio, v 

Radial displacement at r = Rj 

.49 

5.0399xl0- 3 

.499 

5.0602xl0 -3 

.4999 

5.0623xl0 -3 


Figure 8. The thick- walled cylinder problem. Inner radius=3., outer radius=9., 
thickness=l., E— 1000., plane strain condition, mesh=5xl as shown. 
Loading: unit pressure on inner boundary. 

(Consistent units are used for various properties.) 
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Figure 9. The Scordelis-Lo roof problem. Length — 50., radius — 25., thickness — 0.25, 
E=4.32xl0 8 , v=0.3, mesh=NxN. Loading: 90. per unit area in vertical direc- 
tion, i.e.y gravity load; u x =u z = 0 on curved edges. Reference solution: vertical 
displacement at midpoint of free-edge=0.3024 from reference [23]. 

(Consistent units are used for various properties.) 
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